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ABSTRACT 

The existence of planets orbiting a central binary star system immediately raises questions regarding their formation and dynamical 
evolution. Recent discoveries of circumhinary planets by the Kepler space telescope has shown that some of these planets reside close 
to the dynamical stability limit where the strong perturbations induced by the binary makes it very difficult to form planets in situ. 
For binary systems with nearly circular orhits, such as Kepler-38, the observed proximity of planetary orbits to the stability limit can 
be understood by an evolutionary process in which planets form farther out in the calmer environment of the disk and migrate inward 
to their observed position. The Kepler-34 system is different from other systems as it has a high orbital eccentricity of 0.52. Here, we 
analyse evolutionary scenarios for the planet observed around this system using two-dimensional hydrodynamical simulations. The 
highly eccentric binary opens a wide inner hole in the disk which is also eccentric, and displays a slow prograde precession. As a 
result of the large, eccentric inner gap, an embedded planet settles in a final equilibrium position that lies beyond the observed location 
of Kepler-34 b, but has the correct eccentricity. In this configuration the planetary orbit is aligned with the disk in a state of apsidal 
corotation. To account for the closer orbit of Kepler-34 b to the central binary, we considered a two-planet scenario and examined the 
evolution of the system through joint inward migration and capture into mean-motion resonances. When the inner planet orbits inside 
the gap of the disk, planet-planet scattering ensues. While often one object is thrown into a large, highly eccentric orbit, at times the 
system is left with a planet close to the observed orbit, suggesting that Kepler 34 might have had two circumbinary planets where one 
might have been scattered out of the system or into an orbit where it did not transit the central binary during the operation of Kepler. 
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1. Introduction 

With the success of the Kepler space telescope in the past few 
years in detecting several planets around main sequence binaries, 
the formation and dynamical evolution of these objects are now 
among the most outstanding questions in exoplanetary science. 
A survey of the currently known Kepler's circumbinary plane¬ 
tary system^ indicates that, as expected, in all these systems, 
the binaries are close with orbital periods ranging from 7 to 41 
days. The planets in these systems are Saturn-sized or slightly 
smaller and revolve around their host binaries in orbits with pe¬ 
riods from 50 to 300 days. Also, dynamical analyses of these 
objects have indicated that several of these planets are close to 
the boundary of orbital stability (see |Dvorak||1986[ Holman & 
|Wiegert|1999| l, and they are all between two major n : 1 mean- 
motion resonances (MMR). 

The proximity of the orbits of some of the circumbinary 
planets (CBPs) to the stability limit raised the question that 
whether planet formation can proceed efficiently in these re¬ 
gions. Several attempts were made to answer this question 


(Meschiari|2012| Paardekooper et aL|2012| 

Marzari et al. 

2013 

Dunhill & Alexander 20131 |Rafikov||2013 

[Lines et al.[ 

2014 


Ke pler-38 (|Orosz et al.|2012a^ , Kepler-47 i Orosz et al.|2012b i, Kepler- 
64 jSchwamb et al. 2013^ and Kepler-413 ' Kostov et al. 2014 1 . 


* Currently known main sequence binaries with circumbinary planets 
are: Kepler H6 jPoyle et al.|201 f) , Kepler- 34 and 35 ([Welsh et al.|2012[ l. 


not to be possible. The lack of success in growing planetesimals 
to larger sizes in the vicinity of the stability boundary, combined 
with the fact that, as indicated by observations, the orbits of the 
currently known circumbinary planets are almost all coplanar 
with the orbit of their host binaries, strongly suggested that these 
planets formed at large distances (where the perturbing effect of 
the binary on the disk is negligible and planet formation can pro¬ 
ceed in the same way as around single stars) and ended up in 
their current orbits either through planet migration (as a result of 
planet-disk interaction), planet-planet scattering, or a combina¬ 
tion of both ( |Pierens & Nelson|2007||2008b] l. 

The first study of planet migration in circumbinary disks was 
carried out by |Nelson| ( |2003| l who showed that in low eccen¬ 
tricity binaries (ebin ^ 0.2), massive, jovian-type bodies will be 
captured in a 4:1 MMR with the binary and when the binary is 
more eccentric, these planets may be captured in stable orbits 
farther out (see also Pierens & Nelson||2008b| l. In a subsequent 
study, [Pierens & Nelson P007| l extended previous analyses to 
less massive planets (as low as 20 M-Emxh) and showed that mi¬ 
grating circumbinary planets may stop near the edge of the inner 
disk cavity. They suggested that circumbinary planets should be 
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predominantly found in that area. Because in low eccentric bina¬ 
ries, the inner edge of the disk cavity, which is due to the tidal ef¬ 
fect of the binary on the disk material, almost coincides with the 
boundary of orbital stability, the predictions made by these au¬ 
thors agreed with the orbital architecture of the hrst few Kepler 
circumbinary planets. More general cases in which accretion and 
multiple planets were also considered were later studied by the 
same authors ( |Pierens & Nelson|2008a|b] l. 

The earliest attempt to explain the orbital architecture of Ke¬ 
pler CBPs using planet migration was made by |Pierens & Nel- 


son 


(20131. Considering a locally isothermal disk, where the disk 


thermodynamics is modeled by prescribing a given temperature 
profile, and assuming a closed boundary condition at the edge 
of the cavity (i.e., no mass accretion onto the central binary), 
these authors applied their model to the systems of Kepler 16, 34, 
and 35 and showed that the planets in these systems can migrate 
and reside in a stable orbit. However, the values of the semi¬ 
major axes of the final orbits of planets in their models did not 
agree well with their observed values. In a recent paper (Kley & 
Haghighipour 2014), we showed that the limitation of the mod¬ 
els by |Pierens & Nelson] ( |2013| l are due to their two simplifying 
assumptions: considering disk to be locally isothermal and have 
a close boundary condition. A more realistic model requires the 
radiative effects to be taken into account and the disk material to 
flow through the inner edge of the disk into the disk cavity. The 
latter has also been shown in numerical simulations by [Arty 


|mowicz & Lubow| ( |1996 1 and |Gunther & Kl^ ( |2002[ ) where the 


moves the discrepancy observed in the model by Pierens & Nel- 
](2013|l, and allows the planet to reside in close proximity to 


results indicate that despite the appearance of a cavity in the cen¬ 
ter of a circumbinary disk, material can still flow inside and onto 
the central binary. 

In |Kley & Haghighipour ( 2014[ ), we presented an improved 
and extended disk model which included a detailed balance of 
viscous heating and radiative cooling from the surface of the 
disk ( [D’Angelo et al.|2003| l, as well as additional radiative dif¬ 
fusion in the plane of the disk ( |Kley & Crida|200^ . For a planet 
embedded in the disk, this improved thermodynamics can have 
dramatic effects on the planet orbital dynamics such that for low- 
mass planets, it can even reverse the direction of migration ( |Kley| 
& Crida|2008 Kley et al.|2009 1. In addition, we considered an 
open boundary condition and allowed free in-flow of material 
from the disk into the central cavity. This enabled us to construct 
models with net mass in-flow through the disk. We showed that 
when our model is applied to circumbinary planetary systems 
with low eccentricity or circular binaries (e.g., Kepler 38), it re- 


the boundary of orbital stability, between two n : 1 MMRs. 

In our previous study, we considered the system of Kepler- 
38 as this system represents one of the CBP systems with an al¬ 
most circular binary. As mentioned above, in such systems, the 
boundary of the stability of planetary orbits coincides with the 
outer edge of the disk inner cavity. The latter raises the question 
that how the results will change if the binary is eccentric. In a cir¬ 
cumbinary disk around an eccentric binary, the inner cavity will 
be more eccentric. Because the inner cavity is developed due to 
the tidal effect of the binary, its outer edge will no longer co¬ 
incide with the stability boundary as the latter also will change 
with the binary eccentricity. In this paper, we address these is¬ 
sues by applying our disk model to the circumbinary system of 
Kepler 34. The central binary in this system has an eccentricity 
of 0.52 and the orbital eccentricity of its circumbinary planet is 
0.18. It has been suggested that the formation history of this sys¬ 
tem may have been different from the more circular cases such 
as Kepler-16 and Kepler-3 8 ( [Georgakarakos & Eggl|2015] l. 
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This paper is organized as follows: In the next section we 
briefly describe our hydrodynamic model. In section 3, we 
present the structure of the disk without an embedded planet. In 
section 4, the results on the migration of a single planet are pre¬ 
sented which is followed in section 5 by a study on the dynamics 
of a pair of embedded planets. In section 6, we summarize and 
discuss our results. 

2. The hydrodynamic model 

As mentioned in the introduction, the target of our investiga¬ 
tion is the system of Kepler-34, and we taylor our simulations 
to that particular system. Our overall methodology is identical to 
that used in our previous study of the Kepler-38 system ( [Kley &| 
Haghighipour|2014] l, and we give here only a short outline. 

To model the evolution of a disk around a binary star, we as¬ 
sume that the disk is vertically thin and the system is co-planar, 
an assumption well justified by the observations of Kepler cir¬ 
cumbinary planets. We then perform two-dimensional (2D) hy- 
drodynamical simulations in the plane of the binary. To simulate 
the evolution of the disk, the viscous hydrodynamic equations 
are solved in a polar coordinate system (r, 0) with its origin at 
the center of mass of the binary. 

When using the polar coordinate system, there remains a 
hole in the center in which the binary orbits. As pointed out in 
previous studies ( [Pierens & Nelson|2013[pCley & Haghighipouij 
2014| l, the location of the inner boundary of the grid can affect 
the outcome of the simulations. Given that the semi-major axis 
of the central binary of Kepler-34 is ayn = 0.23AU, we consider 
the radial extent of the disk (r) to be from 0.34 AU to 5.0 AU. 
Accordingly, as suggested by studies mentioned above, the in¬ 
ner radius of the grid will be about 1.5 t/bin- The polar angle (p 
varies in an entire annulus of [0,360°]. This domain is covered 
by a grid of 256 x 256 grid cells that are centrally condensed 
in the radial direction and equidistant in azimuth. In all models, 
we evolve the vertically integrated hydrodynamical equations for 
the surface density S and the velocity components (m^, u^). 

When simulating locally isothermal disks, we do not evolve 
the energy equation and instead use an isothermal equation of 
state for the pressure. When considering radiative models, a ver¬ 
tically averaged energy equation is used which evolves the tem¬ 
perature of the midplane ( [Miiller & Kley 2012| l. Radiative ef¬ 
fects are included in two ways. First, a cooling term is con¬ 
sidered to account for the radiative loss from the disk surface 
( [D’Angelo et al.[[2003| l. Second, we include diffusive radiative 
transport in the midplane of the disk using flux-limited diffu¬ 
sion ( [Kley & Crida 2008 [Muller & Kley|201'3] ). In the radiative 


simulations, the ful 
account ([Miiller & 


, vertically averaged dissipation is taken into 
K:iey|2012| l. 


To calculate the necessary height of the disk // at a position 
r, following Gunther & Kley ( 2002|l and Kley & Haghighipour 


P014| l, the individual distances from the two central stars are 
taken into account. The equation of state of the gas in the disk 
is given by the ideal gas law using a mean molecular weight of 
p - 2.35 (in atomic mass units) and an adiabatic exponent of 
j - 1.4. For the shear viscosity we use the a-parametrization 
with a - 0.01 or 0.004 depending on the model, and we set the 
bulk viscosity to zero. For the Rosseland opacity, we use analytic 
formula as provid ed by [Lin & Papaloizou[ ( [1985| l and the flux- 
limiter as in Kley ( 1989^ . 

The calculation of the gravitational force acting by the bi¬ 
nary and planets on the disk is done in the same way as in [Kley 


& Haghighipour (2014jl. Here we use for the planet, a smooth¬ 


ing length of 0.6i/ where the local scale height takes both stars 
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Table 1. The binary parameter and the observed planetary parameter of 
the Kepler-34 system used in the simulations. The mass of the primary 
star is I.OSMq. The values were taken from|Welsh et al.|P012[l. 


a) Binary Parameter 


Mass ratio 

Period 

^bin 

^bin 

q - M 2 /M 1 

[days] 

[AU] 


0.97 

28 

0.23 

0.52 


b) Planet Parameter 


Mass 

Period 

flp 


-^Jup 

[days] 

[AU] 


0.22 

288 

1.09 

0.18 


of the binary into account. In some of the simulations, we add 
a second planet. The integration of the hydrodynamical equa¬ 
tions (explicit/implicit 2nd order scheme) and the N-body inte¬ 
grator (4th order Runge Kutta) are performed again identically 
to Kley & Haghighipour (20141. To calculate the force from the 


disk acting on the planet, we exclude parts of the Hill sphere of 
the planet using a tapering function as given in |Kley & Haghigli^ 
ipour ( |2014| l. We calculate the orbital parameters of the planet 
using Jacobian coordinates, assuming the planet orbits a star 
with mass Mbin - M\ + M 2 , at the binary barycenter. 


2.1. Initial setup and boundary conditions 

In our models, the disk initial surface density is chosen to have a 
S(r) = So profile where r is the distance from the center of 
mass of the binary. For this reference surface density, following 
our earlier results ( |Kley & Haghighipour |2014| ), we chose two 
different values, a higher mass model with So = 3000g/cm^ and a 
lower mass model with a quarter of this value, So/4 = 775g/cm^. 
We chose the initial temperature of the disk to vary with r as 
T(r) cx such that, assuming a central star of Mbin, the ver¬ 
tical thickness of the disk will always maintain the condition 
Hjr - 0.05. For the radiative model, the temperature profile is 
determined self-consistently. The initial angular velocity of the 
disk at a distance r from the barycenter is chosen to be equal to 
the Keplerian velocity at that distance, and the radial velocity is 
set to zero. 

The parameters of the central binary have been adopted from 
Welsh et al. (20121 and are quoted in in Table [T together with 
the planetary data. Owing to the interaction of the binary with 
the disk (and planet), these parameters will slowly vary during 
a simulation. However, during the course of our simulations the 
change is very small. Even for the longest runs that stretch over 
10,000 yrs (equivalent to over 1.3 ■ 10^ binary orbits) the relative 
change of the semi-major axis and the eccentricity of the binary 
is 3 ■ 10“^ and 6 ■ 10“^, respectively. For the models with an 
embedded planet, we consider the planet to have a mass of mp = 
2.1 X 10^'^Mi which is equivalent to mp = 0.22Mjup or about 73 
Earth masses (see Table[T]). At the beginning of each simulation, 
we start the binary at its periastron. 

The boundary conditions of the simulations are constructed 
such that at the outer boundary of the disk, r„ja.x = 5.0 AU, the 
surface density remains constant and equal to its initial value. 
This is achieved by using a damping boundary condition where 
the density is relaxed toward its initial value, l,(r„,ax), and the 
radial velocity is damped toward zero, using the procedure spec¬ 
ified in de Val-Borro et al. ( 2006| l. The angular velocity at the 
outer boundary is also kept equal to its initial Keplerian value. 


Table 2. Parameters of the locally isothermal disk reference model with¬ 
out a planet. 


Parameter 

Value 

Surface density at 1 AU 

3000 g/cc 

Viscosity (a-value) 

0.01 

^min ? ^max 

0.34 AU, 5.0 AU 

fixed H/R 

0.05 



Fig. 1. Azimuthally averaged surface density (top) and disk eccentricity 
(bottom) for our locally isothermal reference model. Shown here are the 
profiles at different evolutionary times (in years); the red curve denotes 
the initial setup. 


Eor the temperature, we use a reflecting condition such that there 
will be no artificial radiative flux through the outer boundary for 
the radiative model. These conditions at the outer boundary lead 
to a disk with zero eccentricity at Hence, r^a.-c has to be 
large enough such that the inner regions are not influenced. 

At the inner boundary (rmm), we consider a boundary con¬ 
dition such that the in-flow of disk material onto the binary is 
allowed. This means, for the radial boundary grid cells at 
we choose a zero-gradient mass out-flow condition, where the 
material can freely leave the grid and flow onto the binary. No 
mass in-flow into a grid is allowed at The zero-gradient con¬ 
dition is also applied to the angular velocity of the material since 
because of the effect of the binary, no well-defined Keplerian 
velocity can be found that could be used otherwise. This zero- 
gradient condition for the angular velocity implies a physically 
more realistic zero-torque boundary. 

With these boundary conditions, the disk can reach a quasi¬ 
stationary state in which there will be a constant mass flow 
through the disk ( |Kley & Haghighipour|2014| l. 
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X 

Fig. 2. Two-dimensional density structure of our isothermal reference 
disk model around the central binary star. The graph shows only a local 
view around the binary. The entire computational grid, however, extends 
from r = 0.34 AU to r = 5.0 AU. The white inner region lies inside the 
computational domain and is not covered by the grid. The positions 
of the primary and secondary stars are indicated by the gray dots. The 
Roche lobe of the secondary is also shown. The central cross marks the 
origin of the coordinate system which coincides with the center of mass 
of the binary. 


3. The structure of the circumbinary disk 

Before we study the evolution of planets within the circumbi¬ 
nary disk around Kepler-34, we analyse the disk’s dynamics due 
to the effect of the central binary. For this purpose, we simulated 
the dynamics of the disk in our reference model (see Table 
with a zero-mass planet. In the simulations the inner boundary 
of the disk is open to mass flow into the central cavity, and at the 
outer boundary the density is relaxed to its initial values. Start¬ 
ing from the initial setup mass flows out into the inner cavity and 
the binary’s torques change the disk’s density profile. However, 
caused by the relaxation outer boundary condition it eventually 
settles into a new final quasi-stationary state, where its mass re¬ 
mains constant and there is a constant mass accretion rate onto 
the central binary. This behaviour is similar to our previous sim¬ 
ulations for the Kepler-38 system ( |Kley & Haghighipour|2014 1 . 
We construct our disk model for a locally isothermal equation 
of state where we do not evolve the energy equation but instead 
leave the temperature of the disk at its initial value. This proce¬ 
dure has the advantage of making the simulations much faster as 
no heating and cooling of the disk have to be considered. 

3.1. Disk structure 

During time, the surface density slowly evolves away from its 
initial profile until it settles in a new equilibrium state. This is 
shown in the upper panel of Fig. where the azimuthally av¬ 
eraged radial surface density profile is shown at different evo¬ 
lutionary times. The equilibration time takes only a few hun¬ 
dred years. After this time, the surface density profiles will no 
longer change considerably. Because of the tidal action of the 
binary on the disk, a central gap is formed with a surface density 
many orders of magnitude lower than inside the disk. The inner 
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Fig. 3. Azimuthally averaged surface density (top) and eccentricity (bot¬ 
tom) for locally isothermal disks, with Hjr = 0.05, around the two sys¬ 
tems, Kepler-34 and Kepler-38. The figure shows the profiles for the 
final equilibrium state without any embedded planet. The radial coor¬ 
dinate has been scaled by the binary separation which is 0.23AU for 
Kepler-34 and 0.15AU for Kepler-38. The results for Kepler-38 have 
been taken from|Kley & Haghighipour]j2014^. 


edge of the disk, which we can define, very approximately, as 
that radius at which the surface density is about half the max¬ 
imum value, lies here at around r 1.25AU. This is about 5 
times larger than abin and as such larger than the stability region 
of massless particle trajectories around binary stars as given by 
|Artymowicz & Lubow| ( T994) l. Our result is in good agreement 
with Pierens & Nelson|(|2013| i and our own results on Kepler-38 
( |Kley & Haghighipour|2014 1 . The larger value of the inner gap 
is caused by the high eccentricity of the binary disk in this case. 
Outside of about r - 3AU, the surface density profile remains at 
the initial value. 


As was shown by [Pierens & Nelson] ( 2013| l, circumbinary 
disks can attain significant eccentricities. We calculate the ec¬ 
centricity of the disk by treating each grid cell as an individual 
particle with a mass and velocity equal to the mass and velocity 
of the cell ( Kley & Dirksen|2006| ). To calculate a radial depen¬ 
dence for the disk eccentricity, edisk(^), we average over the an¬ 
gular direction. In our simulations, the disk eccentricity becomes 
very high, with ejisk about 0.3 to 0.5, in the gap region of the disk, 
as shown in the lower panel of Fig. [T] In the inner, nearly evac¬ 
uated regions, the disk eccentricity can be even higher. Near the 
maximum of the density (at r = 1.7), the eccentricity is around 
0.15 and it drops slowly further out. At radial distances r > 3.0 
AU, the disk eccentricity becomes lower than about 0.01. 


In Fig.[^ we show the 2D surface density distribution for the 
isothermal disk models. We note that this figure shows only the 
inner part of the computational domain around the central binary. 




















































































Kley & Haghighipour: Evolution of circumbinary planets around eccentric binaries: The case of Kepler-34 



30 


25 


-E 20 


15 


10 



0 200 400 600 800 1000 1200 


Time [Years] 


Fig. 4. The evolution of the argument of periapses of the disk with re¬ 
spect to the inertial frame. Results of two different estimates are shown 
for comparison. Red corresponds to estimates using the maximum of 
the 2D density distribution (see Fig.[^, and blue is for estimates using 
the mass weighted mean value of the inner disk. 


Fig. 5. The evolution of the argument of periapses of the binary with 
respect to the inertial frame. The red curve corresponds to the case 
where the backreaction of the disk onto the binary has been articifially 
switched off such that the motion of the binary is not influenced by the 
presence of the disk. The blue curve correponds to the reference model 
as shown in Fig.[^ 


The Roche lobe of the secondary star is shown as well. As shown 
here, an eccentric central binary strongly perturbs the disk and 
produces time varying patterns and a very large eccentric inner 
gap. At the same time, the disk features a strong asymmetric 
maximum in the surface density. 

To test the sensitivity of the results to changes in the inner 
boundary condition, we performed a comparison model with a 
closed inner boundary and found no significant differences in the 
results concerning the gap structure. The width of the gap and the 
disk eccentricity were the same, only the value of the maximum 
density was slightly higher in the closed boundary model. 

In order to study the influence of the binary parameters on the 
disk structure, we compare in Fig. the averaged disk structure 
for the Kepler-34 system and the Kepler-38 system. These two 
systems are different in mass ratio and eccentricity, where for 
the Kepler-38 system, q - 0.26, e = 0.10. It is clear that the 
combination of a higher mass ratio together with a higher binary 
eccentricity creates a much larger variation in the surface density 
profile and a higher disk eccentricity. This difference will have a 
profound effect on the dynamical evolution of embedded planets 
as shown in section |4] 

3.2. Disk precession 

As outlined above, the surface distribution around the central 
binary becomes clearly eccentric. To analyse the disk dynam¬ 
ics without a planet, we further analyzed the precession rate of 
this eccentric mode. In principle, the argument of periapse of the 
disk, mask, can be calculated in the same manner as the disk ec¬ 
centricity, as a mass-weighted average over the disk integrating 
over all annuli ( |Kley & Drrksen|20061 l. However, when calculat¬ 
ing the integral over the entire disk from to r,„ax, the result 
showed a nearly constant value for tiTdisk- However, animations 
of the surface density distribution, as shown for a single snapshot 
in Fig.[^ indicated a clear precession of the inner disk regions. 
For that reason, we decided to restrict the range of integration 
to the inner disk only, and used a radial range extending from 
Tmin = 0.34 to r = 2. As an additional indicator for the preces¬ 
sion, we calculated the disk’s line of periapses from the max¬ 
imum of the density distribution in the disk. The two methods 
resulted in similar results as shown in Fig. I^where the time evo¬ 
lution of Erdi.sk is displayed for the reference model without an 


embedded planet. The above-mentioned averaging method gives 
strong fluctuations whenever the angle crosses zero because in¬ 
dividual rings may have values close to zero or close to 360°. 
The density maximum method on the other hand is noisier due 
to the disk dynamics which results in the ’thicker’ curves. The 
overall agreement of the two methods confirms that the disk ex¬ 
periences a precession. In the initial phase of the simulations 
the precession rate adjusts to the changing density profile and 
in the final equilibrium state the disk experiences a precession 
with a rate of about EJdisk ~ 3°/yr. In comparision to the orbital 
period of the binary, Fbm = 28 d, this indicates a very slow pre¬ 
cession rate. Similar slow precession rates have been found for 
disks around massive planets ( |Kley & Dirksen|2006| ) as well as 
the inner disks in close binary systems such as Cataclysmic bina¬ 
ries ( |Kley et al.|2008| l. The reason that Pierens & Nelson ( |2013| l 
did not find precession of the disk in their models is probably 
related to an integration over the global disk. 


In this paper, we do not investigate further details of our disk 
dynamics, but only briefly comment on the feedback the disk 
can have on the central binary. In Fig.|^ we show the precession 
rate for the binary as induced by the disk. As expected, without 
the disk feedback, the binary orbit does not process. However, 
the inclusion of the disk’s feedack leads to a very slow prograde 
precession of about uTbin ~ 0.05°/yr in the orbit of the binary. 
Here, the total mass of the disk is about 0.015 Mq. During the 
time that the orbit of the binary developed precession, the disk 
yielded a slow decrease of the orbit of the binary with a rate of 
about 4 ■ 10-^AU/yr. 


4. Planet migration 

In this section, we analyse the evolution of embedded planets 
in the disk around Kepler-34. Similar to our previous study on 
Kepler-38 ( |Kley & Haghighipour||2014) l, we first study migra¬ 
tion for our standard isothermal disk model because these are 
numerically faster. Subsequently, we study planets in radiative 
disks, as well. 
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Fig. 6. The evolution of the semi-major axis (top) and eccentricity (hot- 
tom) of an embedded planet in an isothermal disk (reference model) 
that was first brought into an equilibrium. In the two simulations, planet 
was started at different distances from the center of mass of the binary. 
In the first simulation (shown in red), the planet started at a distance of 
flo = 1-75 AU, and in the second simulation (shown in green/blue) it 
started as ao = 2.0 AU. In the simulation shown in red, the planet is im¬ 
mediately released and evolves with the disk. In the second simulation, 
the planet’s orbit is held constant during the first 600 yrs (in blue) and 
then released (green lines). In both simulations, planet migrates inward 
with a rate of about 0.1 AU /lOO yrs. Both models result in unstable 
evolution when the planet has reached a distance of about 1.35 AU. The 
dashed horizontal lines shows the observed semi-major axis and eccen¬ 
tricity of Kepler-34 b. 


4.1. Migration in the standard disk modei 

To study the evolution of the system with an embedded planet, 
we start our simulations with a planet initially placed at dif¬ 
ferent distances (semi-major axes, ao) from the center of mass 
(barycenter) of the binary and in a circular orbit. The planet’s 
evolution is determined by the gravitational action of the disk 
and the central binary. As mentioned earlier, during the evolu¬ 
tion of the planet, its orbital elements are calculated using Ja¬ 
cobian coordinates with respect to the barycenter of the central 
binary. The planet mass is fixed to 73 Earth masses (the observed 
value), and there is no accretion onto the planet. As a reference, 
we present in Table[2 the present orbital parameters of the planet 
Kepler-34 b, as inferred from the observations. 

Figure shows the evolution of the planet through the disk 
for the reference model. We present here the results of two simu¬ 
lations that were carried out for different initial conditions of the 
planet from the center of mass of the binary. In the first model 
(red line), the planet is started directly after insertion. Initially, it 
rapidly migrates inward until it has reached a distance of about 
1.4 AU at around 300 yrs into the evolution. After this, the in- 
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Fig. 7. The evolution of the semi-major axis (top) and eccentricity (bot¬ 
tom) of an embedded planet in isothermal disk models with half the disk 
mass. The red curve refers to the original model and is identical to the 
one shown in Fig.|^ 


ward migration proceeds at a lower pace because the planet has 
reached the eccentric inner hole of the disk. Upon reaching a 
distance of 1.3 AU the eccentricity of the planet has increased to 
about 0.3 and its orbit becomes unstable. To examine whether 
this outcome might have been caused by a too small initial sep¬ 
aration of the planet, so that the system could not reach an equi¬ 
librium, we started a second model with the planet embedded 
initially farther out at ro = 2.0AU. During the first 600 years, 
the semi-major axis of the planet was held constant (blue lines 
in Fig. 1^ which gave enough time to the disk to adjust to the 
presence of the planet. The planet was then released and it mi¬ 
grated due to the effect of the disk torque acting on it (green 
lines). The planet drifted inward with the same speed as in the 
first model, about O.lAU/lOOyrs, and the orbit became unstable 
after the planet reached the same position as before. This implies 
that the outcome of the migration process does not depend on the 
history of the system and is determined solely by the physical 
parameters of the disk. In the following we analyze if reducing 
migration speed can stabilize the orbit. 

4.2. Pianet migration in disks with iower mass 

To investigate whether it was the very rapid inward migration of 
the planet in our standard model that resulted in its unstable evo¬ 
lution, we performed additional isothermal disk simulations with 
a reduced disk mass but otherwise identical parameter. The re¬ 
sults are shown in Fig.[^ Here, the red line refers to the standard 
model as shown with the same color in Fig. and the blue line 
shows a model with 1/2 the disk mass. The reduction of the disk 
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Fig. 8. The evolution of the semi-major axis (top) and eccentricity (bot¬ 
tom) of an embedded planet in disks with lower surface densities (1/4 
of the reference model) and reduced viscosities (a = 0.004). The red 
graph corresponds to the isothermal model and the blue refers to the 
radiative disk. In the two simulations, the planet is started in a circular 
orbit at a distance of Aq = 2.0 AU from the center of mass of the binary. 
In the simulation shown in red, the planet is included in the disk from 
the very beginning of the simulation and it begins to migrate simulta¬ 
neously with the disk evolution starting from its initial density profile. 
In the simulation shown in blue, the disk is first relaxed to the radiative 
equilibrium and then the planet is embedded in it. The dashed horizontal 
lines refer to the observed parameter of the Kepler-34 planet. 


mass resulted in a slower inward migration with slightly reduced 
eccentricity, however, the model became unstable again. 

4.3. Planet migration in disks with lower viscosity and with 
radiation 

In addition to the disk mass, we also varied the viscosity and 
thermodynamics of the disk to investigate its influence on the 
migration process. Fig.j^shows the semi-major axis and eccen¬ 
tricity of the planet embedded in a disk with a lower surface 
density (1/4 of the reference model) and a reduced disk viscos¬ 
ity (a = 0.004). We compare a locally isothermal model (in red) 
with a radiative one (in blue). In the isothermal model, the planet 
migrates inward at a constant rate which is about five times 
slower than in the reference model. When the planet reaches a 
distance slightly inside of r = 1.3 AU, it reverses its direction 
and migrates outward until it finally settles at a distance of about 
1.35 AU. In the radiative case, the planet migrates inward ini¬ 
tially at a fast pace, but then it continues its inward migration 
much more slowly than the isothermal model. The overall evo¬ 
lution of the planet in this case is similar to the isothermal one. 
The final orbit of the planet is only slightly further away at ap¬ 
proximately 1.5 AU. Given the similarity of the isothermal and 



Fig. 9. Graphs of the azimuthally averaged radial surface density of the 
isothermal and radiative disk models of Fig. 8 with an embedded planet. 
The density distributions are taken near the final state of the evolutions 
shown in Fig.[^ The big colored dot in each graph represents the semi¬ 
major axis of the planet at these times. For illustrative purposes, we have 
moved the circles close to their corresponding curves. 


radiative results, we decided to continue our study in the rest 
of the paper by using primarily the isothermal approximation in 
order to cut down on unnecessary computation time. 

Th interesting outward migration of the planet near the end 
of its evolution can be understood in terms of the interaction of 
the eccentric disk with the migrating planet. During the inward 
migration of the planet, the eccentricity of the planet’s orbit re¬ 
mains small. We noticed that when the planet approaches the 
inner, more eccentric regions of the disk, the disk becomes more 
circular which allows further inward migration of the planet. 
Once the planet reaches the inner hole of the disk, the eccentric¬ 
ities of the disk and the planet increase again causing the planet 
move into the disk, periodically. Subsequently, the planet turns 
around and moves slightly outwards, as shown in Fig.[^ 

The final position of the planet in the two models is illus¬ 
trated in Fig. 1^ where the radial surface density distribution is 
plotted together with the location of the planet. In both cases, 
the planet’s final orbit is near the inner edge of the disk where 
the density slope is positive. As shown here, the final position of 
the planet lies just inside of the peak density. This can be com¬ 
pared to the case of Kepler-38 where the central binary is on an 
almost circular orbit and the planet ends up in a position closer 
to the inner binary where the density drops to less than 20% of 
its maximum value ( Kley & Haghighipour|2014| l. 

Fig. 10 shows the two-dimensional density distribution of the 
disk near the final state of planet migration together with the po¬ 
sitions of the binary and planet. Also shown here are the direc¬ 
tions of the periapses of the binary (green arrow), the disk (red), 
and the planet (blue). The blue dashed line refers to the orbit 
of the planet and the white dashed circle to the observed semi¬ 
major axis of the planet in the Kepler-34 system. For the ob¬ 
served orbit, we chose to plot the semi-major axis instead of the 
true ellipse because its orientation is not constant with respect to 
the binary. This figure suggests that the disk’s and the planet’s 
orbits are aligned such that their corresponding periapses always 
point in the same direction. The analysis of the time evolution of 
the system supports this conjecture and confirms that in the equi¬ 
librium state the planet and disk precess exactly with the same 
rate of about 3°/yr (i.e. with the same rate as the disk without 
the planet) see Fig. We also ran a test simulation that contin¬ 
ued this model but with a very small disk mass so that the planet 
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Fig. 10. 2D surface density distribution for the isothermal disk model 
near the final state of its evolution. The gray dots refer to the positions 
of the two binary stars and the planet. The Roche lohe of the system 
of the secondary star and planet are also shown. The blue dashed line 
refers to the orhit of the planet and the blue arrow points to the periapse 
of the planetary orbit. The red arrow points to the periapse of the inner, 
eccentric disk while the green arrow points to the periapse of the binary. 
The white dashed circle refers to the observed semi-major axis of the 
planet with a radius of 1.09 AU. 

would not feel the disk anymore. Interestingly, in this case the 
planet precessed with the same rate as before, indicating that the 
presence of the disk does not have a significant influence on this 
process. This result seems to imply that the precession rate of 
the inner disk equals that of small massless particles, a result 
that will be analysed in our subsequent studies. 

The strict alignment of the planetary orbit with that of the 
eccentric inner disk implies that near apocenter, the planet is al¬ 
ways positioned outside of the maximum density of the disk, as 
indicated with the dotted blue line in Fig.[^ That means, a com¬ 
parison between the semi-major axis of the planet and the radial 
location of the maximum of the azimuthally averaged density 
(see Fig. can give the wrong impression that the planet orbits 
inside the density maximum. 

5. Evolution with two planets 

As shown in previous sections, a single, embedded planet in a 
disk around Kepler 34 stops its migration well outside the ob¬ 
served orbit of Kepler-34 b. This large deviation from the ob¬ 
served orbit is caused mainly by the wide, eccentric inner hole 
of the disk that does not allow the planet to migrate closer to 
the star. One possible solution to this discrepancy would be to 
consider an additional planet in the system located initially at 
a larger distance in the disk. During its inward migration, this 
outer planet, while still embedded in the disk, will interact with 
the inner planet scattering it further inward closer to the central 
binary. 

To examine this scenario, we cari'ied out simulations, start¬ 
ing from the single planet models displayed in Fig. with an 
additional planet with the same mass as that of Kepler-34 b. 
This planet was added to simulations at different evolutionary 
times and was allowed to migrate due to its interaction with the 
disk. All simulations lead to similar results. We show a sample of 
these simulations for an isothermal model in Fig. [TT] In the first 
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Fig. 11. The evolution of the semi-major axis (top) and eccentricity 
(middle) of two embedded planets. The simulation was continued from 
the isothermal model shown in Fig.[^(red line) by adding an additional 
planet with the same mass at 2.0 AU. The green curve corresponds to 
the original model with the new planet added at t x 2300 yrs. The red 
curves correspond to the (original) inner planet and the blue curves is 
for the (additional) outer planet. The bottom panel shows the period 
ratio (outer/inner) of the two planets. 


phase (from t 400 to f2300), the final part of the simulation 
of Fig.j^that was in red, is shown (now in green color) where we 
have shifted the time axis. Then, at f 2300, the second planet 
is added. This new planet migrates inward (blue curve) with an 
initial rate similar to that of the original planet (shown in red). 
At f ss 3200, the migration rate suddenly changes to a slower 
rate. At the same time, the inner planet begins to move inward 
at a similar pace, and its eccentricity is slightly increased (red 
curves). The period ratio Pouter/Anner (lower panel) indicates a 
capture into the 7:4 resonance, which is confirmed by an analy¬ 
ses of the resonant angles. During its subsequent evolution, the 
system remains in resonance and eventually the eccentricity of 
the inner planet reachers high values and the system becomes un¬ 
stable. Similar process in the radiative case of Fig. shows that 
in that simulation, the two planets end up in a 3:2 resonance, and 
again the system becomes unstable. We also carried out simula¬ 
tions where we varied the mass of the additional planet, the disk 
mass, and the time of second planet insertion. We found that in 
general, the system goes through similar evolutionary paths: the 
two planets are captured in different resonances (2:1, 3:2 or 5:3) 
and they ultimately scatter each other in unstable orbits. 
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Fig. 12. The evolution of the semi-major axis (top), eccentricity (mid¬ 
dle) and period ratios of two embedded planets in the disk. The simu¬ 
lation was continued from the isothermal model shown in Fig. (red 
line) that is shown in green here. At t x: 3500yr, a second planet was 
added (blue line) at 2.0 AU. At t x VOOOyr, the planets entered a 3:2 
mean-motion resonance and experienced a scattering event. At a later 
time (f X! 16000), the original evolution with the two planets became 
eventually unstable. At around t x 15000, the disk mass was reduced 
by a factor of 10, and the simulation was continued (shown in purple 
for the outer planet and in light blue for the inner). 


5. 1. Evolution with reduced disk mass 

In the previous section, we noted that the evolution of two plan¬ 
ets embedded in the disk always resulted in their mutual scat¬ 
tering and unstable configurations. Despite this instability, we 
found that after a scattering event, one of the planets remained 
at a location close to the observed orbit of Kepler-34 b. This mo¬ 
tivated us to continue the simulation in one of this cases with a 
reduced disk mass to examine whether the resulted conhguration 
could remain stable. A lower disk mass is to be expected in the 
hnal phase of planet evolution and our mass reduction somewhat 
mimics that state of disk dissipation. 

Fig. [T^ shows the results of one of these simulations. Here, 
we again started the simulation from the isothermal case of 
Fig-Hl now at an earlier time. The evolution of the initial sin¬ 
gle planet is shown again in green. At f « 3500yr into the sim¬ 
ulation, the second planet was added at 2.0 AU (blue line) and 
evolved with the original planet (now shown in red). At around 
t X 7000yr, the system entered a 3:2 MMR followed, shortly, by 
a non-destructive scattering event where the outer planet ended 
up near its starting position at 2 AU and the inner planet moved 


slightly inward, close to the observed location of Kepler-34 b. 
Had we stopped the simulation at this point, we would have had 
a perfect match with the observations for the inner planet. How¬ 
ever, the subsequent evolution (with the disk still being present) 
lead to an inward migration of the outer planet, and a slight out¬ 
ward motion of the inner planet (see appendix |^. At around 
t X 15000, the planets were captured again into the 3:2 reso¬ 
nance which ?At X 16000, resulted in a violent scattering event. 

To study the effect of disk dissipation, we restarted the 
above-mentioned simulation just before the final scattering event 
and after t x 15000, with 1/10 of the original disk mass. For this 
phase, the evolution of the outer planet is shown in light blue and 
that of the inner planet is in purple. The bottom panel shows the 
period ratio in purple as well. During this low disk mass phase, 
the inward migration of the outer planet continues with a very 
small rate, while the inner planet remains at its location. Eventu¬ 
ally the planets are captured in a 3:2 MMR, however, the system 
remains stable as is indicated by the small values of the planets’ 
eccentricities. 


6. Summary and discussion 

Using 2D hydrodynamical simulations, we studied the evolution 
of planets embedded in a circumbinary disk for the parameters of 
the Kepler-34 system. The stellar binary consists of two stars of 
nearly equal masses (~ IMg) orbiting each other in a relatively 
high eccentric orbit with eun - 0.52. These orbital character¬ 
istics of Kepler-34 make this system much more dynamic than 
some of the other systems that host circumbinary planets. 

Our investigation followed that of Kepler-38 (Kley & 
|Haghighipour| |2OT4| | where we adopted a two-step approach. 
We first studied the structure of a circumbinary disk without a 
planet, and then included a planet in the disk and followed its 
subsequent evolution. We considered locally isothermal disks as 
well as more realistic disks that include viscous dissipation, ver¬ 
tical cooling, and radiative diffusion. In the following, we briefly 
present our most important results. 

As shown in Fig. a highly eccentric stellar binary gener¬ 
ates a large inner hole in the disk with a high eccentricity. As a 
result, around an eccentric binary such as Kepler-34, the max¬ 
imum value of the density distribution lies at farther distances 
(around 7 Obm) compared to that in systems with circular binaries 
such Kepler-38 (at about 5 abin)- Our simulations also showed 
that the inner region of the disk will precess in a prograde sense. 
For the system of Kepler-34, the rate of this precession is about 
®7disk = 3°/yr. In turn, the disk that has a total mass of 0.015 
Mq within 5 AU (our reference model), induces a slow prograde 
precession on the binary with a rate of mun - 0.05°/yr. 

To study the evolution of an embedded planet in a disk 
around Kepler-34, we added a planet with a mass of 0.22Mjup 
(the observed value of the mass of Kepler-34 b) to the system 
and allowed the planet to migrate in the disk as a result of disk 
torques acting on it. In our reference disk model, the planet mi¬ 
gration was so rapid that the its orbit became unstable by the time 
it had reached a distance of about 1.35 AU from the barycenter 
of the binary. When reducing the mass of the disk to 25%, and 
considering a lower viscosity (a - 0.004), the orbit of the planet 
became stable. Here, isothermal and radiative models resulted 
in qualitatively similar evolutions. The planet migrated inward 
to a distance of about 1.3 AU from the center and then moved 
out to settle into an equilibrium state of about 1.35 AU for the 
isothermal and 1.5 AU for the radiative models. A test simula¬ 
tion, in which the planet started inside the hole of the disk, at 
the observed distance of 1.09AU, showed that the planet moves 
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outward and settles approximately at the same location as well 
(see Appendix). 

Results of our simulations indicated that the final location of 
the planet is robust and determined by the size and shape (ec¬ 
centricity) of the disk inner hole. The final distance of the planet 
from the binary’s barycenter is larger than the observed semi¬ 
major axis of Kepler-34 b, however its eccentricity is in good 
agreement with its observed value. In this final state, the preces¬ 
sion rate of the planet orbit becomes equal to the disk preces¬ 
sion rate where they enter in a state of apsidal corotation with 
an alignment of their periapses. As a result, the planet always 
orbits outside the disk’s eccentric hole, avoiding moving further 
into the disk. These results seem to indicate that in an eccentric 
binary such as Kepler-34, the evolution of the disk makes it dif¬ 
ficult for the planet to reach close distances from the binary. The 
final position of the planet is determined by the size of the disk 
inner hole which becomes larger for eccentric binaries. 

The difficulty in reaching close orbits for a planet in an ec¬ 
centric disk combined with the discovery of the circumbinary 
planetary system of Kepler-47 with multiple planets motivated 
us to assume that Kepler-34 may contain (or might have con¬ 
tained) at least one additional planet. To investigate the evolution 
of multiple planets in the disk, we added an additional planet, in 
a exterior orbit, to our single-planet disk models and followed 
the combined evolution of both planets. For isothermal and ra¬ 
diative disks, these typically resulted in a capture of the two plan¬ 
ets in low-order, mean-motion resonances with period ratios of 
2:1, 3;2, 5:3 and 7:4. In all cases, the inward migration of the 
resonantly coupled pair resulted in unstable configurations when 
the inner planet reached a distance of about 1.1 to 1.2 AU from 
the binary’s barycenter. In these simulations, one planet was typ¬ 
ically ejected from the system and the second one was scattered 
into an orbit with much larger semi-major axis. 

Upon reducing the disk mass by a factor of 10, we found that 
the planet orbits could stabilize with their subsequent evolution 
resulting in a system where the inner planet has a semimajor axis 
close to its observed value. This suggested that a model in which 
a pair of planets are driven toward the binary just before disk 
dispersal represent a possible scenario for explaining the orbital 
configuration of the Kepler-34 system. In this case, it may be 
interesting to search for additional planets in the system. 

As we have shown in this study, the evolution of planets in 
a disk depends crucially on the structure of the inner hole of the 
disk which itself depends on the binary eccentricity. The detailed 
dependence on disk parameter such as viscosity and disk mass 
will be investigated in more detail in future studies. 

We note that in the models presented in this paper, we re¬ 
stricted ourselves to flat two-dimensional disks. Given the flat¬ 
ness of the observed circumbinary systems this is not an un¬ 
realistic assumption. However, a more comprehensive study of 
the dynamical evolution of these systems requires full three- 
dimensional simulations of radiative disks around binary star 
systems. 
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Appendix A: Evolution of a close in planet 

To test the robustness of the final position of a migrating planet 
in an eccentric disk around the Kepler-34 binary, we performed 
additional simulations where the planet was started inside the 
inner hole of the circumbinary disk. In Fig. |A. 1[ we show the 
evolution of a planet in two systems. In the first system shown 
in red, the planet starts at ro = 2.0 AU as displayed already in 
Fig. In the second case, shown in blue, the planet starts at 
ro = 1.09 AU inside the hole of the disk. In both cases, the initial 
eccentricity of the planet was set to eo = 0. In the second system, 
the starting position of the planet lies at the present observed 
location of Kepler-34 b. 

The simulations show that the planet starting inside the inner 
hole of the disk does not remain at its orbit but moves slowly 
outward. In both systems, the final position of the planet is the 
same, at 1.35 AU from the binary barycenter. The inner planet 
that moves out shows large oscillations in its eccentricity causing 
variations to appear in its semi-major axis. When running simu¬ 
lations for long times, these variations seem to become smaller. 
The outward migration in the second system can be understood 
by the interaction of an eccentric disk with the orbiting planet. 
During one orbit of the planet, because of the oscillations in 
its semi-major axis, it enters periodically into the eccentric disk 
where is experiences a positive net torque that drives the planets 
outward. 



Fig. A.l. The evolution of the semi-major axis of two the planets in an 
isothermal disk. The red curve corresponds to the model displayed in 
Fig. [8] The blue curve corresponds to a planet that has started with an 
initial position of ao = 1.09 inside the inner hole of the disk. 


Article number, page 11 of 1 11 1 









































































